#################################################################
#		Assignment: OneStepWGCNA
#		Author:	Shawn Wang
#		Date: May 14, 2020
#################################################################

##========================Param setting==========================
## Set Param manually
casenum = 20
controlnum = 27
ReadCutoff = 6
SampleNum = 0
title = "zxm384-3w"
RemainGeneNum = 30000

## Set working direction
setwd("/public/home/wangxiao/07.others/ZXM384/02.WorkingDir/")
#################################################################
## Parameters
#' One Step WGCNA
#'
#' @param casenum: case sample number, integer;
#' @param controlnum: control sample number, integer;
#' @param ReadCutoff: Threshold of minimal readcount;
#' @param SampleNum: Threshold of samples number, which samples has the readcount larger than readcount cutoff;
#' @param title: Project Name;
#' @param RemainGeneNum: Retained gene number after the MAD filter step.
#' @notice change the setwd manually.

##====================Section1.Prepration======================
## 02.load package
library(DESeq2)
library(ggplot2)
library(dplyr)
library(WGCNA)
library(stringr)
library(ape)
options(stringsAsFactors = F)
#enableWGCNAThreads() ## for server
## 03.load 
## 3 readcount
load("readcount.Rdata")
## 3 trait
load("trait.Rdata")
## 04.data cleaning
## 4 condition
condition <- factor(c(rep("case",casenum),rep("control",controlnum)),levels = c("case","control"))
colData <- data.frame(row.names = colnames(rawcount), condition)
## remove background noise
x <- rawcount[apply(rawcount,1,function(x) sum(x > ReadCutoff) > (SampleNum*ncol(rawcount))),]
dim(x)
## readcount standardization by DESeq2
dds <- DESeqDataSetFromMatrix(x, colData, design = ~ condition)
dds <- DESeq(dds)
vsd <- assay(varianceStabilizingTransformation(dds))
##====================Section2.WGCNA======================
## 01.parameter
datExpr = data.frame(vsd)
dim(datExpr)
type = "unsigned"
corType = "pearson"
corFnc = ifelse(corType=="pearson", cor, bicor)
maxPOutliers = ifelse(corType=="pearson",1,0.05)
robustY = ifelse(corType=="pearson",T,F)
Title = title
## 1 functions
source("/Users/shawnwang/02.MyScript/OneStepWGCNA/01.Rscript/11.02.WGCNA.SFT.R")
source("/Users/shawnwang/02.MyScript/OneStepWGCNA/01.Rscript/11.03.WGCNA.module.R")
source("/Users/shawnwang/02.MyScript/OneStepWGCNA/01.Rscript/11.04.WGCNA.moduleTrait.R")
source("/Users/shawnwang/02.MyScript/OneStepWGCNA/01.Rscript/11.05.WGCNA.HubGene.R")
## 02. sft
## 2 set remain gene number
GNC = 1-RemainGeneNum/nrow(datExpr) 
## 2 sft calculate
WGCNA.SFT(Title = title,
          datExpr = datExpr,
          GeneNumCut = GNC)
## confirm power
if (is.na(power)){
  power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
                  ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
                         ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
                               ifelse(type == "unsigned", 6, 12))
                 )
  )
}
print(paste("The most suitable Power:",power))

## 2 OneStepNetwork
WGCNA.oneStepNetWork(Title = Title)
## 2 Module-Trait relation
phenotype = pheTmp1
Title = "RawTraitsValue"
WGCNA.ModuleTrait(Title = Title,
                  phenotype = phenotype)
phenotype = pheTmp2
Title = "Tag0-1"
WGCNA.ModuleTrait(Title = Title,
                  phenotype = phenotype)
  
## 2 IntramodularConnectivity
Title = title
connet=abs(cor(datExpr,use="p"))^6
Alldegrees1=intramodularConnectivity(connet, moduleColors)
###(3) Generalizing intramodular connectivity for all genes on the array
datKME=signedKME(datExpr, MEs_col, outputColumnName="MM.")
write.table(datKME, paste(Title,"Conectivity_of_each_modular.xls",sep = "."),
            sep = "\t",
            row.names = T,
            quote = F)
